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ABSTRACT 

Observed molecular clouds often appear to have very low star formation ef- 
ficiencies and lifetimes an order of magnitude longer than their free-fall times. 
Their support is attributed to the random supersonic motions observed in them. 
We study the support of molecular clouds against gravitational collapse by su- 
personic, gas dynamical turbulence using direct numerical simulation. Computa- 
tions with two different algorithms are compared: a particle-based, Lagrangian 
method (SPH), and a grid-based, Eulerian, second-order method (ZEUS). The 
effects of both algorithm and resolution can be studied with this method. We 
find that, under typical molecular cloud conditions, global collapse can indeed 
be prevented, but density enhancements caused by strong shocks nevertheless 
become gravitationally unstable and collapse into dense cores and, presumably, 
stars. The occurance and efficiency of local collapse decreases as the driving 
wave length decreases and the driving strength increases. It appears that local 
collapse can only be prevented entirely with unrealistically short wave length 
driving, but observed core formation rates can be reproduced with more realistic 
driving. At high collapse rates, cores are formed on short time scales in coher- 
ent structures with high efficiency, while at low collapse rates they are scattered 
randomly throughout the region and exhibit considerable age spread. We sug- 
gest that this naturally explains the observed distinction between isolated and 
clustered star formation. 
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1. Motivation 

All presently known star formation occurs in cold molecular clouds. Application of the 
pioneering work of Jeans (1902) on the stability of self-gravitating gaseous systems shows 
that observed molecular clouds vastly exceed the critical mass for gravitational collapse. 
Thus, clouds should efficiently form stars on a free-fall time scale of the order of ~ 10^ 
years in the absence of other effects. However, the lifetime of a typical molecular cloud is 
generally believed to be a factor of ten or twenty longer than predicted by Jeans' classical 
theory (Blitz & Shu 1980), although we note that this is subject to controversy. Ballesteros- 
Paredes, Hartmann, & Vazquez-Semadeni (1999) and Elmegreen (2000) for example have 
argued that not only is the internal structure of molecular clouds transient, but that the 
clouds as a whole may be rather short-lived objects. They suggested that lifetimes of order 
of Tff may actually be necessary to explain the lack of 10 - 20 million year old T Tauri 
stars associated with some molecular clouds. It is observed that stars often do not form 
in one 'catastrophic' event associated with the global collapse of the entire cloud. Instead, 
they form in very localized regions dispersed through an apparently stable cloud (for an 
overview see Williams, Blitz & McKee 1999 and references therein). The total efficiency 
of conversion from gas into stars in typical molecular clouds is very low — of the order of 
a few percent (e.g. Duerr, Imhoff, & Lada 1982; Leisawitz, Bash, & Thaddeus 1989). A 
comprehensive astrophysical explanation remains elusive and its discovery remains one of 
the great challenges for any theory of star formation. 

Molecular clouds are turbulent. This is an essential ingredient for understanding their 
properties and characteristic spatial and temporal behavior. Turbulent gas motions are 
highly supersonic as indicated by the superthermal line widths observed throughout molec- 
ular clouds (Williams et al. 1999). The kinetic energy carried in that motion is sufficient to 
balance the potential energy of the cloud, presumably halting global collapse, a proposition 
that we will test in this paper. However, it can be shown that interstellar turbulence decays 
quite rapidly, on time scales of the order of the free-fall time of the system 



(Mac Low et al. 1998, Stone, Ostriker & Gammie 1998, Mac Low 1999, see also Porter, Pou- 
quet & Woodward 1992a,b, 1994 and Padoan &; Nordlund 1999). Strictly speaking, equation 
|I| is only valid for spherical perturbations with homogeneous density p, with G denoting the 
gravitational constant. However, for more general geometries or density distributions equa- 
tion |I| still gives a good approximation to rfj if we take p to be the mean density of the system. 
To explain their observed long life times, turbulence in molecular clouds must be constantly 
driven (Gammie & Ostriker 1996, Mac Low 1999). The interplay between self-gravity on 
the one hand (leading to local collapse and star formation) and turbulent gas motion on 
the other hand (trying to prevent this process) appears to play a key role in regulating the 
structure of molecular clouds and determining where and when stars form. 
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Unfortunately, a theory of compressible turbulence complete enough to fully address 
the issue of stability against gravitational collapse does not exist, nor is it even visible on 
the horizon. A variety of schemes have been proposed to incorporate the effect of incom- 
pressible turbulence into a gravitational stability analysis. However, molecular clouds are 
extremely compressible. Moreover, the approximations necessary for solution of the resulting 
equations are very stringent and appear to severely limit their applicability to the physical 
conditions found in interstellar clouds. This situation demands a thorough numerical ap- 
proach. Although still far away from fully describing all phenomena present in molecular 
clouds, numerical modeling can capture the important features of supersonic, compressible 
turbulence in self-gravitating, ideal gases. In this paper, we perform a numerical Jeans 
analysis for self-gravitating, compressible, turbulent gas, and apply the result to molecular 
clouds and star forming regions. We do not include magnetic fields here, but the work pre- 
sented here provides the foundation for studies including magnetic fields. Preliminary results 
appear not to reach markedly different conclusions (Mac Low, Klessen, & Heitsch 1999). 

In § ^ we summarize previous work on the question of stability against collapse in self- 
gravitating turbulent media. Then in § ^ we describe our numerical schemes and models. 
The dynamical evolution of our models is discussed in § ^ which also introduces the concept 
of local versus global collapse. In § | we perform a Fourier analysis to quantify the collapse 
behavior on different spatial scales. Section ^ explores the implications of our results for star 
formation in molecular clouds. We speculate about the difference between the 'clustered' 
and 'isolated' modes of star formation and about the different time scales involved. Finally 
in § ^ we summarize our work. 



2. Jeans Analysis 

A first statement about cloud stability can be made from considering the virial theorem. 
Naively speaking, in equilibrium the total kinetic energy in the system adds up to half its 
potential energy, _Ekin + 1/2-Epot = 0. If -Ekin + l/S-Epot < the system collapses, while 
-Ekin + 1/2-Epot > implies expansion. In turbulent clouds, the total kinetic energy includes 
not only the internal energy but also the contributions from turbulent gas motions. If this is 
taken into account, simple energy considerations can already provide a qualitative description 
of the collapse behavior of turbulent self-gravitating media (Bonazzola et al. 1987). 

A more thorough investigation, however, requires a linear stability analysis. For the 
case of an isothermal, infinite, homogeneous, self-gravitating medium at rest (i.e. without 
turbulent motions) Jeans (1902) derived a relation between the oscillation frequency uj and 
the wave number k of small perturbations, 

u;2-c2A;2 + 47rGpo = , (2) 

where Cg is the isothermal sound speed, G the gravitational constant, and po the initial mass 
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density. Note that the derivation includes the ad hoc assumption that the hnearized version 
of the Poisson equation describes only the relation between the perturbed potential and the 
perturbed density, neglecting the potential of the homogeneous solution. This is the so-called 
'Jeans swindle'. The third term in equation ^ is responsible for the existence of decaying 
and growing modes, as pure sound waves stem from the dispersion relation a;^ — Cg/c^ = 0. 
Perturbations are unstable against gravitational contraction if their wave number is below a 
critical value, the Jeans wave number kj, i.e. if 



fc= < fc^^ , lf£A> , (3) 



c: 



or equivalently if the wave length of the perturbation exceeds a critical size given by Aj = 
. This directly translates into a mass limit. All perturbations with masses exceeding 
the Jeans mass, 

.3 _ ^7r\3/2 / 3 



Mj ^ poA^ = [^) Po'ci, (4) 

will collapse under their own weight. As we describe the dynamical evolution of cubic 
subregions inside molecular clouds, we use the cubic definition of the Jeans mass. The 
critical mass for spherical perturbations is lower by a factor of 7r/6. 

Attempts to include the effect of turbulent motions into this analysis were already being 
made in the middle of the century by von Weizsacker (1943, 1951), who also considered the 
production of interstellar clouds from the shocks and density fluctuations in compressible 
turbulence. A more quantitative theory was proposed by Chandrasekhar (1951), who studied 
the effect of microturbulence on gravitational collapse, assuming that collapse occurs on 
scales much greater than the outer scale of turbulence. He derived a dispersion relation 
similar to equation | replacing Cg — > c^+1/3 (f ^) , where (f ^) is the overall velocity dispersion 
due to turbulent motions. Developments through the mid-eighties are reviewed by Scalo 
(1985). Particularly noteworthy is the work of Sasao (1973), who may have been the first 
to quantitatively show that the generation of density enhancements by turbulence, which 
Chandrasekhar (1951) neglected, might be as important as turbulent support. In a more 
recent analysis, Bonazzola et al. (1987) suggested a wave length dependent effective sound 
speed Cg(/c) = + l/3f^(/c), leading to a dispersion relation 

cu^ - (cl + Iv^k)^ e + 47rG'po = . (5) 

In this description, the stability of the system depends not only on the total amount of 
energy, but also on the wave length distribution of the energy, since v^{k) depends on the 
turbulent power spectrum e{k) as 

POO 

v\k) = / e{k')dk' . (6) 



k 



Thus, the system can be stable at some wave lengths, but not at others. This approach was 
also adopted by Vazquez- Semadeni & Gazol (1995), who added Larson's (1981) empirical 
scaling relations to the analysis. 
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The most elaborate investigation of the stability of turbulent, self- gravitating gas was 
made by Bonazzola et al. (1992), who used renormalization group theory to derive a dis- 
persion relation with a generalized, wave number- dependent, effective sound speed and an 
effective kinetic viscosity that together account for turbulence at all wave lengths shorter 
than the one in question. They found a general dispersion relation (their equation 4.13) that, 
if applied to turbulence with a power-law energy spectrum e{k) = Ak'", predicts a critical 
value of the power-law exponent a = 3. According to their analysis, turbulence with a spec- 
trum steeper than this can support a region against collapse at large scales, and below the 
thermal Jeans scale, but not in between. On the other hand, they claim that turbulence with 
a shallower slope, as is expected for incompressible turbulence (Kolmogorov 1941), Burgers 
turbulence (Lesieur 1997), or shock dominated flows (Passot, Pouquet & Woodward 1988), 
cannot support clouds against collapse at scales larger than the thermal Jeans wave length. 

These analytical approaches make a strong assumption that substantially limits their 
reliability, namely that the equilibrium state is homogeneous, with constant density po- How- 
ever, observations clearly show that molecular clouds are extremely non-uniform. In fact, it 
may even be possible to describe the equilibrium state as an inherently inhomogeneous ther- 
modynamic critical point (de Vega, Sanchez & Combes 1996a,b; de Vega & Sanchez 1999). 
As a consequence of the assumption of homogeneity, the further assumption of microturbu- 
lence must then be made. The largest turbulent scale is significantly smaller than the scale 
of the analysis. Interstellar turbulence, however, does not appear to exhibit such a cut-off in 
the power spectrum, but rather extends over all spatial scales present in the system. A fur- 
ther corollary of the assumption of homogeneity is that the turbulent dynamical time scale 
is much shorter than the collapse time scale Tg, which is only justified if the assumption of 
microturbulence holds. 

One way to achieve progress and circumvent the restrictions of a purely analytical ap- 
proach is to perform numerical simulations. Bonazzola et al. (1987), for example, used low 
resolution (32 x 32 collocation points) calculations with a 2-dimensional spectral code to 
support their analytical results. Also restricted to two dimensions were the hydrodynami- 
cal studies by Passot et al. (1988), Lcorat, Passot & Pouquet (1990), Vazquez-Semadeni et 
al. (1995) and Ballesteros-Parcdcs, Vazquez-Semadeni & Scalo (1999), although performed 
with far higher resolution. Magnetic fields were introduced in two dimensions by Passot, 
Vazquez-Semadeni, & Pouquet (1995), and extended to three dimensions with self-gravity 
(though at only 64^ resolution) by Vazquez-Semadeni, Passot, & Pouquet (1996). A careful 
analysis of 1-dimensional computations including both MHD and self-gravity was presented 
by Gammie & Ostriker (1996), who extended their work to 2.5 dimensions more recently 
(Ostriker, Gammie, & Stone 1999). Preliminary results of high-resolution (256^ zone) sim- 
ulations with MHD and self-gravity have been presented by Mac Low et al. (1999) and by 
Ostriker, Gammie, & Stone (1998). In the present paper we use two numerical algorithms to 
examine the stability properties of three-dimensional hydrodynamical turbulence at higher 
resolution than before. In subsequent work in preparation we will include magnetic fields as 



- 7- 



well. 



3. Numerical Methods 

Direct numerical simulation of the Euler equations for gas flow does not reach the 
enormous Reynolds numbers typical of molecular clouds due to the intrinsic or numerical 
viscosity of any finite-difference or spectral method. However, if the details of behavior at 
the dissipation scale do not affect the behavior of larger scales, then all that is required is a 
low enough viscosity to separate the two scales. Incompressible turbulence appears to behave 
like this (e.g. Lesieur 1997). Resolution studies of energy decay in supersonic compressible 
turbulence suggest that it may also be true in this case (Mac Low et al. 1998a). The 
resolution studies we do here also address this question, as increasing the resolution decreases 
the dissipation scale, which is always close to the zone size. 

We use both Lagrangian and Eulcrian numerical methods to solve the equations of 
self-gravitating hydrodynamics in three dimensions (3D) in an attempt to bracket reality 
by taking advantage of the strengths of each approach. This also gives us some protection 
against interpreting numerical artifacts as physical effects. The Lagrangian method we use 
is smoothed particle hydrodynamics (SPH), while the Eulerian method is the astrophysical 
hydrocode ZEUS. In future work we use this numerical calibration in the interpretation of 
self- gravitating MHD models computed with ZEUS. 

3.1. SPH 

SPH is a Lagrangian, particle-based scheme to solve the equations of hydrodynamics. 
The fluid is represented by an ensemble of particles, each carrying mass, momentum, and 
thermodynamical properties. Fluid properties at any point are obtained by averaging over 
a set of neighboring particles. The time evolution of the fluid is represented by the time 
evolution of the particles, governed by the equation of motion and equations to implement 
hydrodynamic properties. The technique can therefore be seen as an extension of the pure 
gravitational A^-body system. Excellent overviews of the method, its numerical implemen- 
tation, and some of its applications are given by the reviews by Benz (1990) and Monaghan 
(1992). The code used here derives from a version originally developed by Benz (1990). 
It includes a standard von Neumann-type artificial viscosity (Monaghan & Gingold 1983) 
with the parameters — 1 and (5y — 2 for the linear and quadratic terms. The system is 
integrated in time using a second-order Runge-Kutta-Fehlberg scheme, allowing individual 
time steps for each particle. Furthermore, the smoothing volume over which hydrodynamic 
quantities are averaged in the code is freely adjustable in space and time such that the 
number of neighbors for each particle remains approximately fifty. 
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SPH can resolve very high density contrasts because it increases the particle concen- 
tration, and thus the effective spatial resolution, in regions of high density, making it well 
suited for computing collapse problems. Conversely, it resolves low-density regions poorly. 
Shock structures tend to be broadened by the averaging kernel in the absence of adaptive 
techniques. It is also very difficult to include magnetic fields in the algorithm. SPH can 
be run on the special-purpose hardware device GRAPE (Sugimoto et al. 1990, Ebisuzaki et 
al. 1993; and also Umemura et al. 1993, Steinmetz 1996), which allows supercomputer-level 
calculations to be done on a normal workstation. As we concentrate on subregions inside 
molecular clouds of much larger extent, we use periodic boundary conditions, as implemented 
by Klessen (1997) on GRAPE. 

The correct numerical treatment of gravitational collapse requires the resolution of the 
local Jeans mass at every stage of the collapse (Bate & Burkert 1997). In the current code, 
once an object with density beyond the resolution limit of the code has formed in the center 
of a collapsing gas clump it is replaced by a 'sink' particle (Bate, Bonnell, & Price 1995). This 
particle has a fixed radius on the order of the Jeans length at the threshold density. We set 
this density to be 10^ times the average density in the simulation, which roughly corresonds 
to the maximum resolvable density contrast. The sink particle inherits the combined mass 
of the replaced SPH particles, as well as their linear and angular momenta. It has the ability 
to accrete further SPH particles from its infalling gaseous envelope, which are then removed 
from the computation. Adequately replacing high-density cores and keeping track of their 
further evolution in a consistent way prevents the time step from becoming prohibitively 
small. We are thus able to follow the collapse of a large number of cores until the overall gas 
reservoir becomes exhausted. 



3.2. ZEUS-3D 

ZEUS-3D is a well-tested, Eulerian, finite-difference code (Stone & Norman 1992, Clarke 
1994). It uses second-order van Leer (1977) advection, and resolves shocks using von Neu- 
mann artificial viscosity. Self-gravity is implemented via an FFT-solver for Cartesian coordi- 
nates. It also includes magnetic fields in the magnetohydrodynamic approximation. For the 
models discussed here, we use a three-dimensional, periodic, uniform, Cartesian grid. This 
gives us equal resolution in all regions, and allows us to resolve shocks well everywhere. On 
the other hand, collapsing regions cannot be followed to scales less than one or two cells. 

We must again consider the resolution required for gravitational collapse. For a grid- 
based simulation, the criterion given by Truelove et al. (1997) holds. Equivalent to the 
SPH resolution criterion, the mass contained in one grid zone has to be smaller than the 
local Jeans mass throughout the computation. Applying this criterion strictly would limit 
our simulations to the very first stages of collapse, as we have not implemented anything 
like sink particles in ZEUS. We have therefore extended our models beyond the point of 
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full resolution of the collapse, as we are primarily interested in the formation of collapsed 
regions, but not their subsequent evolution. Thus, in the ZEUS models, the fixed spatial 
resolution of the grid implies that strongly collapsed cores have a larger cross-section than 
appropriate for their mass. In encounters with shock fronts the probability for these cores 
to get destroyed or lose material is overestimated. Cores simulated with ZEUS are therefore 
more easily disrupted than they would be physically. SPH, on the other hand, underestimates 
the disruption probability, because sink particles cannot lose mass or dissolve again once they 
have formed. The physical result is thus bracketed by these two numerical methods (also see 



3.3. Models 

We perform our computations using normalized units. The considered volumes are cubes 
with side L = 2, extending from -1 to 1, which are subject to periodic boundary conditions 
in every direction. The total mass in the box is M = 1, therefore the uniform initial density 
is po = 1/8- We use an isothermal equation of state, with sound speed Cg = 0.1, chosen to 
set the number of thermal Jeans masses contained in the box to A^j = 64. Time is measured 
in units of the initial global free-fall time of the system. 

To generate and maintain turbulent flows we introduce Gaussian velocity fluctuations 
with power only in a narrow interval k — 1 < \k\ < k, where k = L/Xd counts the number 
of driving wave lengths in the box. This offers a simple approximation to driving by 
mechanisms that act on that scale. Comparing runs with different k will then give some 
information on how, for example, turbulence driven by large-scale shearing motions might 
differ from turbulence driven by low-mass protostars. We set up the initial velocity field as 
described in Mac Low et al. (1998), with perturbations drawn from a Gaussian random field 
determined by its power distribution in Fourier space, following the usual procedure. For 
each three-dimensional wave number k we randomly select an amplitude from a Gaussian 
distribution around unity and a phase between zero and 2tt. We then transform the resulting 
field back into real space to get a velocity component in each zone, and multiply by the 
amplitude required to get the desired initial root mean square (rms) velocity on the grid. 
We repeat this for each velocity component independently to get the full velocity field. 

To drive the turbulence, we then normalize this fixed pattern to produce a set of pertur- 
bations (5z7(x, y, z), and at every time step add a velocity field 6v{x, y, z) = ASuto the velocity 
V. The amplitude A is chosen to maintain constant kinetic energy input rate Ei^ = AE/At. 
For a compressible flow with a time- dependent density distribution, we maintain a constant 
energy input rate by solving a quadratic equation in the amplitude A at each time step, as 
discussed in Mac Low (1999). In dynamical equilibrium, the driving luminosity Ei^ equals 
the rate of turbulent energy dissipation. To estimate the input rate necessary to reach and 
maintain a certain equilibrium level of the kinetic energy we use equation 7 of Mac Low 
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(1999). We find that this equation underestimates the driving energy needed to maintain 
the SPH models at a specific equilibrium kinetic energy by 20-30% for reasons that we do 
not yet fully understand. Comparisons with other techniques will probably be required to 
resolve this discrepancy. We drive the SPH models somewhat harder to compensate, as can 
be seen m table m by comparing Ei^ for SPH and ZEUS models with equivalent driving wave 
length fcdrv and E^f^. Dynamical equilibrium is reached typically within one global shock 
crossing time t = L/{v). The equilibrium value is determined to an accuracy of better than 
10%. Keeping the energy input unaltered we then switch on self-gravity with gravitational 
constant G = 1, and allow the evolution to proceed. This defines t = in our models. Their 
most important properties are summarized in table 0. 



3.4. Scaling 

The dynamical behavior of isothermal self-gravitating gas is scale free and depends only 
on the ratio between potential energy and kinetic energy (including thermal energy). We 
can scale our models to physical units with a mass scale of the thermal Jeans mass Mj given 
by equation (H), a length scale given by the Jeans length Aj derived from equation d), and 
a time scale given by the free-fall time scaled from equation (|I|) 

-»=(0-34Myr)(^5jiL_)"^ (7) 

where G = 6.67 x 10~® cm^ g~^ s~^, and the number density is taken as = p/fi, with 
yU = 2.36 mn. The 120 thermal Jeans masses in our simulation cubes then correspond to 

M = (413M0) ( ' ( , \ , (8) 

^ °^ VO.2 km s"V V104 cm-3y ' ^ ^ 

where a sound speed Cg = 0.2 km s~^ corresponds to a temperature T = 11.4 K with the 
value of fi we use. Finally we may compute the size of our cube by noting that the Jeans 
length in our computational units is Aj = O.ly^Svr ^ 0.501, while the size of the cube is 
L = 2, so that in physical units 

L = (0.89 pc) ( ( / A . (9) 

^ ^ ' VO.2 km s"V \W cm-y ^ ' 

As an example, let us consider a dark cloud like Taurus with n{E.2) ~ 10^ cm"^, and 
Co ~ 0.2 km s^^. Then our simulation cube holds a mass M = 4.1 X WMq and has a size 
L = 8.9 pc. The time unit (free-fall time scale) is tq = 3.4 Myr, and the average thermal Jeans 
mass for the homogeneous distribution follows as Mj = 65 Mq . Another example would be 
a dense cloud forming massive stars such as the BN region in Orion, with n{B.2) ~ 10^ cm"^ 
and Cs ~ 0.2 kms~^. Here the simulated cube holds a mass of M = 130 and is of size 
L = 0.28 pc. The time unit is now = 0.1 Myr, and the thermal Jeans mass is Mj = 2.1 Mq. 
(Note again, that in the spherical definition the Jeans mass is smaller by a factor 7r/6.) 
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4. Local vs. Global Collapse 

In this section we begin by showing numerical results that suggest that local collapse 
can occur in turbulent self-gravitating media even if the kinetic energy contained in the 
system is sufficient to stabilize it on global scales (§ The strong shocks ubiquitous in 

supersonic turbulence compress small regions sufficiently that the turbulence can no longer 
support them. We then consider what promotes or prevents this process (§ ^3), and the 
importance of turbulent collapse in real molecular clouds (§ |4.3| ). 

4.1. Local Collapse in a Globally Stable Region 

We compute models with both SPH and ZEUS in which the turbulence is driven at 
strengths above and below the critical value needed to prevent gravitational collapse accord- 
ing to the analytic predictions of § ^. The models can be characterized by two parameters, 
the kinetic energy before gravity is turned on, and the typical driving wave number k at 
which energy is injected (see § p.3[ ). We define an effective turbulent Jeans mass (Mj)turb by 
substituting Cg — > + 1/3 (f^) for the thermal sound speed Cs in equation where we 
approximate the rms velocity of the flow {v^) by 2£'kin/M. We do simulations with (Mj)turb 
of 0.6, 3.2, and 18.2. These values have to be compared to the total system mass M = 1 
in order to determine whether global stability is reached. Note that our deflnition of the 
Jeans mass uses the mean density in the simulations. This is equivalent to examining the 
collapse properties of isolated gas cubes. In inflnite media (local) density contrasts should 
be used instead. For this reason, the quoted turbulent Jeans masses are lower limits to the 
true critical values for support against gravitational collapse. The true stabilizing effect of 
turbulence on large scales is stronger than indicated from merely comparing (Mj)turb with 
the total mass in the system. 

We flnd that local collapse occurs even when the turbulent velocity fleld carries enough 
energy to counterbalance gravitational contraction on global scales. This confirms the re- 
sults of two-dimensional (2D) and low-resolution (64^) 3D computations with and without 
magnetic fields by Vazquez- Semadeni et al. (1996). An example of local collapse in a globally 
supported cloud is given in figure |I[ It shows a sequence of 3D density cubes of the SPH 
model B2h which is driven in the wave length interval 3 < /c < 4 so that the turbulent Jeans 
mass (Mj)turb = 3.2. The first cube shows the system at t = 0.0. Hydrodynamic turbulence 
is fully established but gravity has not yet been included in the computation. (Note again 
that time is measured in units of the global free-fall time of the system rg and the zero-point 
is set when gravity is switched on). The second cube shows the system at a time t = 1.1. 
Density fiuctuations generated by supersonic turbulence in converging and interacting shock 
fronts that locally exceed the Jeans limit begin to contract. The central regions of these 
high-density clumps have undergone sufficient gravitational contraction to be identified as 
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collapsed cores. Numerically, in the SPH code they have been replaced with sink particles. 
There are altogether twelve cores containing = 5% of the total gas mass in the system. 
At t = 3.9 the number of dense embedded cores has grown to 46 and they account for 25% 
of the mass. At t = 7.1 roughly 50% of the gas mass is accreted onto 53 dense cores. The 
first cores form in small groups randomly dispersed throughout the volume. Their velocities 
directly reflect the turbulent velocity field of the gas they are created from, in which they 
are still embedded, and from which they continue to accrete. However, as more and more 
mass accumulates on the cores the gravitational interaction between the cores themselves 
increasingly determines their dynamical evolution. The core cluster begins to behave more 
like a collisional A^-body system, in which close encounters are dynamically important. 

Local collapse in a globally stabilized cloud is not predicted by the analytic models 
described in Sec. |^. For the parameters of the models presented here, the dispersion rela- 
tion equation ^ forbids gravitational contraction at any scale. However, this equation was 
derived under the assumption of incompressibility. The presence of shocks in supersonic 
turbulence drastically alters the result, as was first noted by Elmegreen (1993) and studied 
numerically by Vazquez-Semadeni et al. (1996). The density contrast in isothermal shocks 
scales quadratically with the Mach number, so the shocks driven by supersonic turbulence 
create density enhancements with 5p o: M.^, where M. is the rms Mach number of the flow. 
In such fluctuations the local Jeans mass is decreased by a factor of and therefore the 
likelihood for gravitational collapse increased. 

To test this explanation numerically, we designed a test case driven at short enough 
wave length and high enough power to support even fluctuations with 5p (x ^A'^, and ran it 
with both codes. To ensure sufficient numerical resolution for these models, and P5, we 
computed a subvolume of mass M = 0.25 with reduced sound speed Cg = 0.05 driven at wave 
number /c = 9 — 10. This is equivalent to an effective driving wave number = 39 — 40 on the 
regular cube (M = 1, Cg = 0.1). Within 20 rg neither of these models show signs of collapse. 
All the other globally supported models with less extreme parameters that we computed did 
form dense cores during the course of their evolution, supporting our hypothesis that local 
collapse is caused by the density fluctuations resulting from supersonic turbulence. 

The two numerical methods that we use are complementary, as discussed in § ^. SPH is 
a particle based, Lagrangian scheme. It resolves regions of high density well, and the use of 
sink particles makes it straightforward to define dense cores, but it does not resolve shocks 
well. Once a collapsing region passes beyond the density threshold and is converted into a 
sink particle, it cannot be destroyed. It continues to accrete matter from its surroundings 
and to interact gravitationally with other cores. This overestimates the survival probability 
of collapsing, Jeans- unstable fiuctuations. ZEUS, on the other hand, is an Eulerian grid 
method, well suited for resolving shocks, but worse at modeling gravitational collapse. Due 
to the fixed grid, it overestimates the volume of collapsed cores, leading to an enhanced 
cross-section to destructive processes such as tidal interactions between cores, or the pertur- 
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bations of passing shock fronts. Hence, the probabihty for core formation and survival in 
the turbulent environment is underestimated. The real behavior of self-gravitating, turbulent 
gas lies in between, bracketed by the two methods which we apply here. Both methods show 
local collapse occurring in globally stabilized clouds. 

Figure ^ illustrates this point by comparing 2D slices through 3D models which are run 
by the two different codes with similar turbulent driving power spectra at both medium and 
high resolution. Each slice is centered on the densest core on the grid. (Because we use 
periodic boundary conditions, we are free to shift the window across the simulated volume 
in any direction. These boundaries do not introduce artificial perturbations.) Figures |^a and 
b show the SPH models B2 and B2h with 50 000 and 200 000 particles, while figures ^ and d 
show the ZEUS models P2 and V2h with 128^ and 256^ grid zones. We use a new realization 
of the initial conditions with the same statistical properties for each of these models, so there 
is no expectation that they will have identical structures, only that they will have similar 
typical structures. The roundish appearance of structures in the SPH models, especially at 
lower resolution, stems from the smoothing intrinsic to the SPH algorithm. The Lagrangian 
nature of the scheme leads to high spatial resolution in high-density regions but degraded 
resolution in low-density regions where particles are sparse. Conversely, ZEUS does very 
well at modeling the shock and void structure, especially in the high-resolution model V2h, 
but the dense collapsed cores are underresolved. The shocks and filaments clearly resolved 
by the ZEUS model are also present in the SPH model, but tend to be rather smeared out 
by the lack of resolution in the lower density regions. Nevertheless, all the images clearly 
indicate the presence of strong shocks which sweep up gas into gravitationally collapsing 
regions. 



4.2. Promotion and Prevention of Local Collapse 

The total mass and lifetime of a fluctuation determine whether it will actually collapse. 
Roughly speaking, the lifetime of a clump is determined by the interval between two succes- 
sive passing shocks: the first creates it, while if the second is strong enough, it disrupts the 
clump again if it has not already collapsed (Klein, McKee & Colella 1994, Mac Low et al. 
1994). If its lifetime is long enough, a Jeans unstable clump can contract to sufficiently high 
densities to effectively decouple from the ambient gas flow. It then becomes able to survive 
the encounter with further shock fronts (e.g. Krebs & Hillebrandt 1983), and continues to 
accrete from the surrounding gas, forming a dense core. The weaker the passing shocks, and 
the greater the separation between them, the more likely that collapse will occur. Equiva- 
lently, weak driving and long typical driving wave lengths enhance collapse. The influence 
of the driving wave length is enhanced because individual shocks sweep up more mass when 
the typical wave length is longer, so density enhancements resulting from the interaction of 
shocked layers will have larger masses, and so are more likely to exceed their local Jeans 
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limit. Turbulent driving mechanisms that act on large scales will produce large coherent 
structures (filaments of compressed gas with embedded dense cores) on relatively short time 
scales compared to small-scale driving even if the total kinetic energy in the system is the 
same. 

We demonstrate the effect of the driving wave length in figure which compares SPH 
model Blh with driving wave numbers = 1 — 2 to model B3 driven with = 7 — 8 at 
a time when sink particles have accreted 5% of the gas mass. (These density cubes can 
be directly compared with figure |l]b, which shows the intermediate case = 3 — 4 at the 
same evolutionary stage.) Note the difference in the morphology of the density structures. 
Figure ^ is dominated by one large shock front that traverses the volume, which is the sole 
site of core formation. On the other hand, the density structure in model B3 (figure ^) 
is far more homogeneous, without any large-scale structure. Cores form alone, randomly 
dispersed throughout the volume. This comparison is discussed below in § |6.1| . 

The influence of driving strength and wave length on local collapse can be examined 
by measuring the amount of mass accreted onto collapsed regions over time in each model. 
In the SPH models, this can be computed quite simply by adding up the masses of the 
sink particles at each time. As ZEUS does not include sink particles, we instead employ 
a modified version of the Clumpfind method (Williams, de Geus, & Blitz 1994; see also 
appendix 1 in Klessen & Burkert 2000). In this routine, clumps are defined as regions of 
connected zones whose densities lie above a certain threshold. In order to be able to use 
Clumpfind on models as large as 256^ zones, we replaced the inefficient clump identification 
routines with an algorithm based on the dilation operators implemented in IDL. We use two 
criteria to separate collapsed cores from shock-generated fluctuations. First, we require the 
average density of those cores to exceed the mean value expected for isothermal shocks, 
p > Ai^po- Here, Ai is the rms Mach number and po is the mean density. Second, we count 
only fluctuations for which the potential energy exceeds the kinetic energy, E^°^'^ < \E^'^^\, 
and which are more massive than the local Jeans mass, M^, > Mj°^^(p). We use logarithmic 
density contours instead of linear ones in order to get a wide enough density range so that 
most detected clumps consist of more than one cell. However, we also accept single high- 
density cells as cores. These are common at late stages of the evolution, when the envelopes 
of cores have been removed by further shock interactions and only the collapsed centers 
remain. In this case the we use the ratio of potential to internal energy as the criterion for 
collapse. 

A figure of merit that we can use to examine the effect of driving strength and wave 
length is the time ^5% needed to sweep up 5% of the mass into compact cores. Table |l| 
describes several sequences of models identical except for their driving wave length: the high 
and medium resolution SPH models Al — A3, Bl — B4 and C2, and the low, medium, and 
high resolution versions of the ZEUS models VI — V3. Comparison of the values of for 
these models shows that collapse and accretion occurs more rapidly for models with larger 
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driving wave length (smaller driving wave numbers and larger typical scales). Comparison 
of the SPH models with k = 3 — 4 shows that stronger driving also delays collapse. 

Models run with the same driving strength at different resolutions and with the two 
different codes can be compared to determine the level of numerical convergence, and the 
effect of the different algorithms. Comparison of SPH models B2£ to B2h shows that a change 
of linear resolution of 2.2 yields a change in t^x of only 12.5%. Similarly comparison of the 
128^ zone to the 256^ zone resolution ZEUS models shows better than 10% agreement, except 
for high wave number driving, where the disagreement is still less than 25%. Comparison of 
the ZEUS and SPH models with k = 3—4 driving of the same strength also shows better than 
25% quantitative agreement. The 2D cuts through medium and high resolution models with 
both codes with the same driving wave length {k = 3 — 4) and driving strength ((Mj)turb = 
3.2) presented in figure ^ visually demonstrate the level of morphological agreement. We 
emphasize that the qualitative result that local collapse occurs at a rate dependent on the 
driving wave length and strength is recovered at all resolutions and with both codes. 

A more detailed understanding of how local collapse proceeds comes from examining 
the full time history of accretion for each model. Figure ^ shows the accretion history for 
three sets of SPH models. For each set of models, the driving strength is held constant while 
the effective driving wave length is varied, showing the pronounced effect of the wave length 
at equal driving strength. At the extreme, if the driving is at wave lengths below the Jeans 
wave length of the shocked layers local collapse does not occur (model i35). The A models 
have lower driving strength than the B and C models, demonstrating the effect of driving 
strength at each driving wave length. 

The cessation of strong accretion onto cores occurs long before all gas has been accreted. 
This appears to be because the time that dense cores spend in shock-compressed, high-density 
regions decreases with increasing driving wave number and increasing driving strength. In 
the case of long wave length driving, cores form coherently in high-density regions associated 
with one or two large shock fronts that can accumulate a considerable fraction of the total 
mass of the system. The overall accretion rate is high and cores spend suffient time in this 
environment to accrete a large fraction of the total mass in the region. Any further mass 
growth has to occur from chance encounters with other dense regions. In the case of short 
wave length driving, the network of shocks is tightly knit. Cores form in shock generated 
clumps of small masses because individual shocks are not able to sweep up much matter. 
Furthermore, in this rapidely changing environment the time interval between the formation 
of clumps and their destruction is short. The period during which individual cores are 
located in high-density regions where they are able to accrete at high rate is short as well. 
Altogether, the global accretion rates are small and begin to saturate at lower values of 
as the driving wave length is decreased. 

Figure ^ shows the accretion history for the three 256^ ZEUS models, Vlh to 'D3h. The 
fractional core mass in the model with large scale driving {T>lh) is strongly affected by 
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the large shocks that run through the volume. At t ^ 2.1, for example, a shock destroys 
the most massive core, so drops suddenly. Between successive shock passages, the cores 
recover, so they accrete a substantial mass fraction over the run. Models V2h and V3h with 
k = 3 — 4 and k = 7 — 8 display a steady mass growth similar to the SPH models. The 
more frequent shocks in these models reduce the accretion rate by stripping away material 
from the vicinity of the central high-density zones. These isolated zones do not lose mass 
from shock encounters, but are subject to numerical clipping, so the measured fraction is, 
as explained before, a lower limit to the actual accretion fraction. A clear indication of local 
collapse is once again seen. 

To further investigate the influence of numerical resolution, figure ^ compares the time 
history of accretion for SPH models with varying particle numbers, but identical turbulent 
Jeans mass (Mj)turb = 3.2 and driving wave numbers 3 < A; < 4 (see table |I|). The difference 
in effective linear resolution (cube root of the particle number) between B2i and B2h is 2.2. 
We also had to distinguish the effects of statistical variance from the effects of resolution. 
To do this, we repeated the intermediate resolution simulation B2 four more times, varying 
only the random seeds used to generate the Gaussian fields (models B2°--B2'^, dashed lines). 
We actually find stronger variation between the different models at the same resolutions 
than between models at different resolutions, suggesting that numerical diffusivity doesn't 
have as large an effect as the natural statistical variation. This is not surprising given the 
stochastic nature of turbulent flows. Protostellar cores form in molecular clouds through 
a sequence of highly probabilistic events. Especially at late times, their mass accretion 
is strongly influenced by chaotic A^-body dynamics (Klessen et al. 1998, Klessen Sz Burkert 
2000). All models agree well at early times when initial local collapse occurs, suggesting that 
we are well converged on our basic result. At late times, variations between the different 
models become stronger. These variance effects need to be kept in mind when interpreting 
the accretion rates of individual models. For the ensemble average at late times we do not 
expect significant variations at the different numerical resolutions which we study, though 
our current set of calculations is not large enough to quantify this statement. 

Figure ^ shows the time history of accretion for ZEUS models with the same parameters 
as the SPH models shown in figure ^, and numerical resolution increasing from 64'^ to 256'^ 
zones. Strong fluctuations in the lower resolution curves are caused by core disruptions due 
to shocks, which cannot occur in the SPH models as sink particles are never destroyed. The 
fluctuations decrease with increasing resolution because the cores have smaller cross section 
in the high resolution models, and are thus less liable to be destroyed by shocks. The high 
resolution model T>2h shows a well defined accretion behavior and reaches a saturation level 
at a mass fraction of about 8%, where the local free-fall time of the cores is roughly equal 
to the time interval between two shock passages. All three models reach this level at least 
intermittently, suggesting it defines a reasonably firm upper limit for these ZEUS models, 
and thus a lower limit to the amount of mass that can actually be accreted under these 
physical conditions, with the SPH models giving an upper limit. 
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4.3. Application to Molecular Clouds 

The global star formation efficiency in normal molecular clouds is usually estimated be 
of the order of a few per cent. Their life times are typically thought to be a fewxlO^ years 
which is equivalent to a few tens of their free-fall time tq (Blitz & Shu 1980, Blitz 1993, 
Williams et al. 1999). It would be consistent with these estimations if the mass fraction 
of protostellar cores in our simulations remained below M^, = 5% for 10 rg. Indeed, as 
indicated in column 8 of table |l| local collapse can be slowed down considerably in the 
case of small-scale driving. However, if the hypothesis of rapid molecular cloud evolution 
is correct (Ballesteros-Paredes et al. 1999, Elmegreen 2000), the constraints on the driving 
scale and strength are substantially changed. Furthermore, it needs to be noted that the 
local star formation efficiency in molecular clouds can reach very high values. For example, 
the Trapezium star cluster in Orion is likely to have formed with an efficiency of about 
50% (Hillenbrand & Hartmann 1998). In § |6.1| we will argue that this apparent difference 
between the 'clustered' and 'isolated' model of star formation can be explained in terms of 
the properties of the underlying turbulent velocity field of the parental gas. 

The energy dissipation scale in molecular clouds should also be considered. It was ffist 
shown by Zweibel & Josafatsson (1983) that ambipolar diffusion would be the most important 
dissipation mechanism in typical molecular clouds with very low ionization fractions x = 
Pi/Pn, where pi is the density of ions, p„ is the density of neutrals, and p = pi + pn- An 
ambipolar diffusion strength can be defined as 

>^AD = v\/Uni, (10) 

where v\ = B^/inpn approximates the effective Alfven speed for the coupled neutrals and 
ions if Pn ^ Pi, and u^i = 'ypi is the rate at which each neutral is hit by ions. The coupling 
constant depends on the cross-section for ion-neutral interaction, and for typical molecular 
cloud conditions has a value of 7 ~ 9.2 x 10^^ cm^ s~^ g~^ (e.g. Smith & Mac Low 1997). 
Zweibel & Brandenburg (1997) define an ambipolar diffusion Reynolds number as 

Rad = LV/\ad = MALUni/vA, (11) 

which must fall below unity for ambipolar diffusion to be important, where L and V are 
the characteristic length and velocity scales, and /Aa = V/va is the characteristic Alfven 
Mach number. In our situation we again can take the rms velocity as typical value for V. 
By setting Rad = 1, we can derive a critical length scale below which ambipolar diffusion is 
important 

Lcr = -r^ ^ (0.041 pc) 3 ,^ , (12) 

where the magnetic field strength B = 10-Bio pG, the ionization fraction x = 10~^X6, the 
neutral number density = lO^ris cm~^, and we have taken p„ = /in„, with p = 2.36 mu. 
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We can attempt to compare this value to the numerical dissipation scale by directly com- 
puting the ratio of the thermal Jeans length Aj that we use to scale our models (as discussed 
in § ^.41 ) to Lcr- We do this by assuming that ionization and magnetic field both depend on 
the density of the region, following the empirical laws rij = 3 x 10~^ cm~^ (r;,„/10^ cm~^)^/^ 

1 /2 

(e.g. Mouschovias 1991), and Biq ~ ?>n.^ (e.g. the observational summary of Crutcher 
1999). We can then find the interesting result that 

T^ = 16.1-MAf-— (13) 



Lcr \0.2 km s 

Crutcher (1999) suggested that typical values of the Alfven Mach number M.a are only 
slightly above unity. With the value in our simulations Aj = 0.1-\/87r 0.5 (cf. §3.4) and 
noting that our cube has a side length of L = 2 ^ 4Aj, this implies that the critical length 
scale on which ambipolar diffusion becomes important in our model units is = L/64. 
This is comparable to or even slightly greater than the length on which numerical dissipation 
acts in our highest resolution models. Thus, we can conclude somewhat surprisingly that we 
may be close to capturing the full dissipation-free range available to real molecular clouds 
in our models. 



5. Fourier Analysis 

In this section we discuss the energy distribution on different spatial scales during various 
stages of the dynamical evolution of the system. We perform a Fourier analysis of the 
energy, computing the power spectra of kinetic and potential energies. To allow for a direct 
comparison, all models are analyzed on a Cartesian grid with 128^ cells. For the SPH 
models this is done using the kernel smoothing algorithm, while the 256^-ZEUS models 
are simply degraded in resolution. For each cell the potential and kinetic energy content 
is calculated, and the kinetic energy is further decomposed into its solenoidal (rotational) 
and compressional parts. These quantities are then all transformed into Fourier space, to 
find the contribution of different dimensionless wave numbers k, or equivalently, to find the 
distribution of energy over different spatial scales A^ = L/k. 

The energy spectrum of fully developed turbulence for small-, medium- and large-scale 
driving is shown in figure |. It shows the SPH models (a) Al, (b) A2 and (c) ^3 just at 
the time t = 0.0 when gravity is turned on. In each plot the thick solid lines describe the 
potential energy as a function of wave number k, and the thick long-dashed lines represent the 
kinetic energy, which can be decomposed into its solenoidal (rotational) and compressional 
components. They are defined via the velocities by V -Vsoi = and V x -Ucom = 0, respectively. 

The models Al and A2, which are driven at long and intermediate wave lengths {k = 
1 — 2 and k = 3 — 4), appear to exhibit an inertial range below the driving scale, i.e. between 
0.5 ^ logiQ k ^ 1.5. Note that, in real clouds, the dissipation scale may lie near the upper 



- 19 - 



end of this wave number range as discussed in § O. In this interval the energy distribution 
approximately follows a power law very similar to that predicted by the Kolmogorov (1941) 
theory (-Ekin k~^/^). This is understandable given that, in our models, once turbulence is 
fully established, the solenoidal component of the kinetic energy always dominates over the 
compressible one, Eso\ > Ecom- For a pure shock dominated flow (-Ecom ^ -£^soi) one would 
expect a power spectrum with slope —2 (Passot et al. 1988). To guide the eye, both slopes 
are indicated as thin dotted lines in plots (a) to (c). For model ^3 the smaller number of 
available modes between the driving scale k = 7 — 8 and the Nyquist frequency does not 
allow for an unambiguous identification of a turbulent inertial range. The permanent energy 
input necessary to sustain an equilibrium state of turbulence produces a signature in the 
energy distribution at the driving wave length. This is most clearly visible in figure 

The system is globally stable against gravitational collapse, as indicated by the fact that 
for every wave number k the kinetic energy exceeds the potential energy. For comparison we 
plot in figure ||d the energy distribution of a system without turbulent support. The data 
are taken from Klessen et al. (1998) and stem from an SPH simulation with 500 000 particles 
containing 220 thermal Jeans masses and no turbulent velocity field, but otherwise identical 
physical parameters. The snapshot is taken at t = 0.2rff after the start of the simulation. 
This system contracts on all scales and forms stars at very high rate within a few free-fall 
times Tfj. Contrary to the case of hydrodynamic turbulence, the kinetic energy distribution 
is dominated by compressional modes, especially at small wave numbers. The overall energy 
budget is determined by the potential energy -Epot, which outweighs the kinetic energy E^^i^ 
on all spatial scales k by about an order of magnitude. 

Figure || concentrates on model B2h with (Mj)turb = 3.2 and = 3 — 4. It describes 
the time evolution of the energy distribution. Figure ^ shows the state of fully established 
turbulence for this model just when gravity is turned on (t = 0.0). In the subsequent 
evolution, local collapse occurs in shock-generated density enhancements where the potential 
energy dominates over the kinetic energy. This affects the small scales first, as seen in the 
plotted time sequence. As collapse progresses to higher and higher densities, the scale 
where the potential energy dominates rapidly grows. Once the mass fraction in dense cores 
has reached about ~ 3%, the potential energy outweighs the kinetic energy on all scales. 
However, this should not be confused with the signature of global collapse. The power 
spectrum of the potential energy is constant for all k. It is the Fourier transform of a delta 
function. Local collapse has produced point-like high-density cores. The overall filling factor 
of collapsing clumps and cores is very low, so most of the volume is dominated by essentially 
pure hydrodynamic turbulence. As a consequence, the velocity field on large scales is not 
modified much (besides a shift to higher energies). On small scales, however, the flow is 
strongly influenced by the presence of collapsed cores which is noticeable as a flattening 
of the power spectra at large wave numbers. Despite their small volume filling factor, the 
cores dominate the overall power spectrum. The solenoidal part of the kinetic energy always 
dominates over the compressional modes and also the signature of the driving source in the 
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energy spectrum remains, visible as a 'bump' in the kinetic energy spectrum at A; ~ 8. 

To show that the global features of our models are insensitive to the numerical method 



used, in figure |10| we compare the energy spectra of four different simulations with identical 
physical parameters. Analogously to figure H, we chose simulations B2, B2h, T>2 and V2h 
which all have (Mj)turb = 3.2 and A; = 3 — 4. Models B2 and B2h are SPH simulations 
with 50 000 and 200 000 particles, while V2 and T>2h were calculated using the ZEUS code 
with a resolution of 128^ and 256^ grid zones, respectively. Figures |TOa-d directly compare 



the different energy components in the four models at t = 0, at the stage of fully developed 
pure hydrodynamic turbulence just before gravity is switched on. The sequence |10|e-h does 
the same after gravity has been switched on and the first collapsed cores have formed at 
^5%, when the mass accumulated in dense cores is M^, = 5% of the total mass. This state is 
identical to the one depicted in figure 0, allowing for direct comparison. 

Comparing the spectra of the different models during the stage of pure hydrodynamic 
turbulence (figures p!0|a-d) shows excellent agreement between the energy spectra of the 
different models, suggesting the energy distribution is well converged. Between the scales 
of energy input (at = 3 — 4) and diffusive energy loss, all the spectra follow the same 
power law with slope —5/3 (analogous to the spectra shown in figure |^). The dissipation 
scale manifests itself as a drop-off from the power-law at large wave numbers. The inertial 
range of turbulence is largest in the high-resolution ZEUS model V2h, where it spans about 
one order of magnitude in k. The high-resolution SPH and medium-resolution ZEUS models 
B2h and T>2 have inertial ranges nearly as long. The medium-resolution SPH model B2 has 
the shortest range with Alog^g^ ~ O-^- wave numbers above the dissipation scale, our 
results are converged to better than 10% in the log of energy. 

In the presence of self-gravity, the energy spectra are no longer well converged. The 
actual density contrast reachable in collapsing cores, or to a lesser extent in shock fronts, 
depends on the numerical resolution and algorithm used (see § |). The same applies to the 
potential energy and to the compressional component of the kinetic energy. Figures |TO| e and 
h therefore exhibit significant differences between the various models. These differences are 
much smaller for the solenoidal part of the kinetic energy, which measures rotational motions 
and is therefore less sensitive to strong density contrasts in small volumes. Variations in the 
total kinetic energy distribution are mainly due to differences in the compressional modes. 

The rapid energy decrease for wave numbers log^^Q A; > 1.4 in the grid-based model V2 
is due to the fact that these scales approach the grid resolution. A similar decrease would 
be seen in the other three models if they were sampled at wave numbers all the way up 
to the effective resolution (grid size for ZEUS or smoothing length for SPH). Remember 
that all spectra shown here are computed on the same grid with a linear resolution of 128 
cells. Despite the fact that the 256^ ZEUS model V2h has been resampled and degraded 
in resolution, large density contrasts still occur on the smallest scales of the resampled grid. 
The energy spectra therefore remain fiat towards the Nyquist wave number. Similarly, the 
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use of adaptive particle smoothing lengths in SPH allows the resolution of dense cores smaller 
than the cell size of the 128^ grid used for the energy sampling. Again, there is no loss of 
power towards the Nyquist wave number of the spectra. However, high resolution in high 
density regions is achieved at the cost of low resolution in voids. As low density regions 
occupy most of the volume, on large scales the SPH simulations tend to have lower energy 
content than the grid-based models. 

Understanding these various effects allows us to understand what questions we can 
reasonably ask of these simulations. The presence or absence of collapse and the distribution 
of kinetic energy on large scales are questions for which we can give well-converged answers, 
but the details of the strength of that collapse still depend on the details of the numerical 
method and should not be used quantitatively. 

6. Implications for Star Formation in Molecular Clouds 

In § ^ we have shown that the rate and efficiency of local collapse in turbulent molecular 
clouds depend on the strength and the effective wave length of the driving energy input. Star 
formation will follow local collapse (e.g. Vazquez-Semadeni, Canto, & Lizano 1998), so we 
can use these properties of our turbulence models to try to explain the observed spatial and 
age distributions of young stars in molecular clouds. We use the spatial and age distributions 
of sink particles generated in the SPH models with different parameters for this purpose. 

6.1. Clustered vs. Isolated Star Formation 

Different star formation regions present different distributions of protostars and pre- 
main sequence stars. In some regions, such as the Taurus molecular cloud, stars form isolated 
from other stars, scattered throughout the cloud (Mizuno et al. 1995). In other regions, they 
form in clusters, as in L1630 in Orion (Lada 1992), or even more extremely in starburst 
regions such as 30 Doradus (Walborn et al. 1999). 

/^From the simulations presented here, it is evident that the length scale and strength at 
which energy is inserted into the system determine the structure of the turbulent flow and 
therefore the locations at which stars are most likely to form. Large-scale driving leads to 
large coherent shock structures (see e.g. flgure |^). Local collapse occurs predominantly in 
fllaments and layers of shocked gas and is very efficient in converting gas into stars. This 
leads to what we can identify as 'clustered' mode of star formation: stars form in coherent 
aggregates and clusters. Even more so, this applies to regions of molecular gas that have 
become decoupled from energy input. As turbulence decays, these regions begin to contract 
and form dense clusters of stars with very high efficiency on about a free-fall time scale 
(Klessen et al. 1998, Klessen & Burkert 2000). The same holds for insufficient support. 
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i.e. for regions where energy input is not strong enough to completely balance gravity. They 
too will contract to form dense stellar clusters. 

The 'isolated' mode of star formation occurs in molecular cloud regions that are sup- 
ported by driving sources that act on small scales and in an incoherent or stochastic manner. 
In this case, individual shock induced density fluctuations form at random locations and 
evolve more or less independently of each other. The resulting stellar population is widely 
dispersed throughout the cloud and, as collapsing clumps are exposed to frequent shock 
interaction, the overall star formation rate is low. 

To demonstrate these points, we compare in figure |TT] the distribution of sink particles 
for several different models, projected onto the xy- and xz-planes. As an example of co- 
herent local collapse, we choose model Bl, where the turbulence is driven strongly at long 
wave lengths. The flow is dominated by large coherent shocks, so cores form in aggregates 
associated with the filamentary structure of shock compressed gas (cf. with figure 0). The 
overall efficiency of converting gas into stars in this 'clustered' mode is very high. The upper 
half of figure |11] compares the model Bl with model B3, which is driven at small scales and 
results in incoherent collapse behavior. Individual cores form independently of each other 
at random locations and random times. In this 'isolated' mode, cores are widely distributed 
throughout the entire volume and exhibit considerable age spread. 



In the lower half of figure |TT| we contrast the large-scale driving model Bl with a simu- 
lation of freely decaying turbulence described by Klessen (2000) that has the same thermal 
Jeans mass. In decaying turbulence, once the kinetic energy level has decreased sufficiently, 
all spatial modes of the system contract gravitationally, including the global ones. As in 
the case of large-scale shock compression, stars form more or less coevally in a very limited 
volume with high efficiency. Both insufficient turbulent support and the complete loss of it 
therefore appear to lead to clustered star formation. The Trapezium cluster in Orion may 
be a good example for the outcome of this mechanism (e.g. Hillenbrand 1997, Hillenbrand 
& Hartmann 1998). All the projections shown in figure ^ are taken at a stage of the dy- 
namical evolution when the mass accumulated in dense cores is M^, ^ 20%. This occurs at 
very different times, as noted in the captions, which directly reflects the varying efficiencies 
of local collapse in these models. 

Despite the fact that turbulence that is driven on large scales as well as turbulence 



that is freely decaying lead to star formation in aggregates and clusters, flgure suggests a 
possible way to distinguish between them by taking the long-term evolution of the resulting 
clusters into account. Whereas decaying turbulence typically leads to the formation of a 
bound stellar cluster, the dynamical relaxation of aggregates associated with large-scale 
coherent shock fronts quite likely results in their complete dispersal. This is illustrated in 



flgure 0, which compares the core distribution in model Bl and in the decay simulation at 
M* fa 65%, when both systems have already undergone considerable evolution. The cores in 
model Bl are completely dispersed throughout the molecular cloud volume. The cluster that 



- 23 - 



formed during the turbulent decay remains bound with a much longer evaporation time scale. 
Note, however, that at late stages of the dynamical evolution our isothermal model becomes 
less appropriate as the feedback effects from newly formed stars are not taken into account. 
Ionization and outflows from these stars will likely retard or even prevent the accretion of 
the remaining gas onto the protostars. This limits the validity of our models at very late 
times. 



6.2. The Time Scales of Star Formation 

In the previous section we have argued that stellar clusters form predominantly in 
molecular cloud regions that are insufficiently supported by turbulence or where only large- 
scale driving is active. In the absence of driving, molecular cloud turbulence decays more 
quickly than the free-fall time scale (Mac Low 1999). The free-fall time tq is thus the typical 
time scale on which dense stellar clusters will form in the absence of support (Klessen et 
al. 1998) or in the presence of decaying turbulence (Klessen 2000). Even in the presence of 
support from large-scale driving, collapse will occur on roughly this time scale, as shown for 



model Bl in figure |T3|a. If we assume that once we have identified a dense core it continues 
to collapse on a very short time scale to build up a stellar object in its center, then this 
spread relates directly to the star formation time scale. Therefore the age distribution will 
be roughly Tg for stellar clusters that form coherently with high star formation efficiency. 
When scaled to low densities (n(H2) ~ lO^cm"^ and T ^ lOK) the global free-fall time 
scale in our models is rfj = 3.3 x 10^ years. If star forming clouds such as Taurus indeed 
have ages of order rg, as suggested by Ballesteros-Paredes et al. (1999), then the long star 
formation time scales computed here is quite consistent with the very low star formation 
efficiencies seen in Taurus (e.g. Leisawitz et al. 1989), as the cloud simply has not had time 
to form many stars. In the case of high- density regions (n(H2) ~ 10^ cm~'^ and T 10 K) the 
dynamical evolution proceeds much faster and the corresponding free-fall time scale drops 
to Tfj = 1.0 X 10^ years. These values indeed agree well with observational data, e.g. the 
formation time scale of the Orion Trapezium cluster, which is inferred to stem from gas of 
density n(H2) ^ lO^cm"^, is estimated to be less than 10^ years (Hillenbrand & Hartmann 
1998). 

The age spread increases with increasing driving wave number k and increasing (Mj)turb- 
Molecular cloud regions supported against global collapse by driving sources that act on small 
scales host stochastical star formation on much longer time scales and with much lower 
efficiency. The process is incoherent and the expected stellar age spread therefore larger. 
Indeed, in figure |l3|b, which shows the accretion history of selected cores in model B2 with 
k = 3 — A (representing more isolated star formation), core formation extends over a longer 
period. This is even more pronounced in model B3 with k = 7 — 8 shown in figure p!3|c. Note 
that the real time spread in this model is even larger than suggested by the figure, because 
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by the time we stopped the simulation the accreted mass fraction was only = 35%. We 
expect that more cores would form in the subsequent evolution. Models Bl and B2, on the 
other hand, already reach M^, ^ 70% in the time interval shown. They each form roughly 50 
cores, twice as much as model B3. For a direct comparison, figure |TB|d plots the distribution 
of core formation times in each of the three models on the same scale. These long periods 
of core formation for globally supported clouds appear consistent with the low efficiencies of 
star-formation in regions of isolated star formation, such as Taurus, even if they are rather 
young objects with ages of order rg. 

7. Summary and Conclusions 

We have studied the conditions that allow self gravity to cause collapse in a region of 
supersonic turbulence. We used this study to determine whether interstellar turbulence can 
support molecular clouds against gravitational collapse, reveahng the scales and physical 
conditions that allow star formation to occur. To perform these studies, we computed 
numerical simulations of the time evolution of turbulent, self-gravitating, isothermal gas 
with two different computational schemes: a particle-based, Lagrangian method (SPH); and 
a second-order, Eulerian, grid-based method (ZEUS). By comparing results from these two 
different numerical schemes we benefit from the advantages of both methods, and we are 
furthermore able to estimate the infiuence of algorithm as well as resolution on our results. 
We next summarize and discuss our results. 

1. Supersonic turbulence strong enough to globally support a molecular cloud against 
collapse will usually cause local collapse. The turbulence establishes a complex net- 
work of interacting shocks. The local density enhancements in fiuctuations created by 
converging shock fiows can be large enough to become gravitationally unstable and 
collapse. This occurs if the local Jeans length becomes smaller than the size of the 
fiuctuation. The probability for this to happen, the efficiency of the process, and the 
rate of continuing accretion onto collapsed cores are strongly dependent on the driving 
wave length and on the rms velocity of the turbulent flow, and thus on the driving 
luminosity. Collapse criteria derived from incompressible, self-gravitating turbulence 
(Chandrasekhar 1951, Bonazzola et al. 1987, 1992, Vazquez- Semadeni & Gazol 1995) 
indeed determine the global or large-scale collapse properties of the medium. However, 
the occurrence and ubiquity of local collapse in shock-generated fluctuations drastically 
limit the application of these criteria to interpreting the actual behavior of star-forming 
regions, as localized collapse can still occur even if the cloud as a whole is stabilized by 
turbulence. 

2. Fluctuations in turbulent velocity flelds are highly transient. The random flow that 
creates local density enhancements can also disperse them. For local collapse to result 
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in the formation of stars, locally Jeans unstable, shock-generated, density fluctuations 
must collapse to sufficiently high densities on time scales shorter than the typical time 
interval between two successive shock passages. Only then are they able to 'decouple' 
from the ambient flow pattern and survive subsequent shock interactions. (If they begin 
collapse magnetically supercritical, they will remain so for the rest of the collapse.) 
The shorter the time between shock passages, the less likely these fluctuations are 
to survive. Hence, keeping the scale of energy input fixed and increasing the driving 
luminosity leads to a decrease of the star formation efficiency. Local collapse takes 
longer to occur and the mass accretion rate onto cores is reduced. Similarly, driving 
on small scales leads to a lower star formation rate than driving on larger scales at the 
same rms velocity. Quantitatively, our models appear to show that it is possible to 
prevent 95% of the gas from collapsing into dense cores over ten global free-fall times 
with strong enough driving on short enough wave lengths. If a physical mechanism 
for such driving can be found, this could indeed explain the long cloud life times and 
low star formation rates commonly ascribed to Galactic molecular clouds (Blitz & Shu 
1980, Bhtz 1993). Conversely, if such driving does not exist, then molecular clouds 
should be transient objects and the short life times proposed by Ballesteros-Paredes et 
al. (1999) and Elmegreen (2000) appear more likely. 

3. Local collapse can only be halted completely if the turbulent driving mechanism sup- 
plies enough energy on scales smaller than the Jeans length of the 'typical' fiuctuation. 
In supersonic turbulence the typical density contrast is 5p oc A4^, where A4 is the rms 
Mach number of the fiow. Thus, the Jeans length is reduced by a factor of A4 with re- 
spect to the global value. Complete prevention of local collapse requires even stronger 
and shorter wave length driving, as there will be stochastic turbulent fiuctuations with 
even larger density contrast. However, the time scale for the occurence of high density 
fluctuations increases rapidly with 6p, so sufliciently strong driving can prevent local 
collapse for arbitrarily long periods of time. Such strong driving may be rather difficult 
to arrange in a real molecular cloud, however. 

If we assume that stellar driving sources have an effective wave length close to their 
separation, then the condition that driving acts on scales smaller then the Jeans wave 
length in 'typical' shock generated gas clumps requires the presence of an extraor- 
dinarily large number of stars evenly distributed throughout the cloud, with typical 
separation 0.1 pc in Taurus, or only 350 AU in Orion (taking our fully supported case 
as an example). This is not observed. Mac Low et al. (1999) show that magnetic 
fields probably cannot transfer energy efficiently enough to small scales either. Fur- 
thermore, ambipolar diffusion may begin to damp turbulent motions at these scales. 
Unless some other mechanism can force energy onto these small scales, local collapse 
will occur within globally supported molecular clouds. In addition, very small driving 
scales seem to be at odds with the observed velocity fields at least in some molecular 
clouds (e.g. Heyer et al. 1997 for the Cep 0B3 cloud). 
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4. Interstellar clouds driven on large scales or without even global turbulent support very 
rapidly form stars in clusters. Gas collapses into dense cores within a few free-fall 
times and the star formation efficiency is more than 50%. On the contrary, in gas that 
is supported by turbulence driven at small scales, local collapse occurs sporadically 
over a large time interval, forming isolated stars. The total star formation efficiency 
before the cloud dissolves due to stellar feedback or external shocks will probably be 
low. Thus, the strength and nature of the turbulence may be fully sufficient to explain 
the difference between the observed isolated and clustered modes of star formation. 

5. In turbulent flows, it is impossible to predict from the start when and where individual 
cores form and how they evolve. Firm statistical results can, however, be derived from 
analyzing large ensembles of cores and from characterizing other global indicators of 
the dynamical state of the system such as total potential and kinetic energy. In all 
our models except the ones driven below the fluctuation Jeans scale, gravity eventually 
begins to dominate over kinetic energy. This flrst occurs on small scales, indicating 
the presence of local collapse. As dense collapsed cores form, the power spectrum of 
the gravitational energy becomes essentially flat. The kinetic energy, on the other 
hand, appears to follow at intermediate wave numbers a Kolmogorov power spectrum 
with slope —5/3, less steep than the spectrum expected for pure shock flows. The 
slope remains almost unaltered during the course of the evolution, indicating that 
a large volume fraction of the system is always well described by pure hydrodynamic 
turbulence. The spatial extent of collapsing regions (where infall motions dominate over 
the turbulent flow) is relatively small. This also explains the fact that the solenoidal 
component of the flow always dominates over the compressional part. 
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Table 1: Overview of the models. The columns give model name, numerical method, 
resolution, driving wave lengths k, energy input rate Ei^^, equilibrium value of kinetic energy 
without self-gravity E^^^, turbulent Jeans mass (Mj)turbj and the time required to reach a 
core mass fraction = 5%. The resolution is given for SPH as particle number and for 
ZEUS as number of grid cells. Dashes in the last column indicate that no sign of local 
collapse was observed within 20Tff, while stars indicate that the numerical resolution was 
insufficient for unambiguous identification of collapsed cores. The total mass in the system 
is M = 1. Models B5 and V5 focus on a subvolume with mass M — 0.25 and decreased 
sound speed Cg = 0.05. They are driven with A; = 9 — 10 and E'in = 0.06. When scaled up to 
the standard cube this corresponds to the effective values given in square brackets. Model B2 
has been calculated five times with different random initializations. The additional models 
are not listed separately, but are called B2"' - B2'^ in the text. 





Fig. 1. — SPH density cubes for model B2h, which is driven in the interval 3 < A; < 4, 
shown (a) at the time when gravity is turned on, (b) when the first dense cores are formed 
and have accreted M* = 5% of the mass, (c) when the mass in dense cores is M* = 25%, and 
(d) when = 50%. Time is measured in units of the global system free-fall time scale Tg. 
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Fig. 2. — Comparison of 2D density slices through 3D models with identical physical pa- 
rameters ((Mj)turb = 3.2 and A; = 3 — 4) computed with different numerical methods and 
resolution: SPH models (a) B2 and (b) B2h with 50 000 and 200 000 particles, and ZEUS 
models (c) V2 and (d) V2h with 128^ and 256'^ grid cells. For further details see table |I[ 
To allow for comparison, the time is chosen such that the mass accreted onto dense cores 
is = 5%. Density is scaled logarithmically with the separation of contour levels being 
one decade. Each cut is centered on the density maximum in the simulation. In SPH, the 
density distribution has been interpolated onto a uniform grid using kernel smoothing. The 
arrows indicate the velocity comoonents in the olane of section. 
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Fig. 3. — Density cubes for models (a) Blh {k = 1 — 2) and (b) B3 {k = 7 — 8) at dynamical 
stages where the core mass fraction is M^, = 5%. Compare these figures with figure [l|b. 
Together they show the influence of different driving wave lengths for otherwise identical 
physical parameters. Note the different visual appearance of the systems and the different 
times at which M^, = 5% is reached. 
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Fig. 4. — Fraction of mass in dense cores as function of time. All models are computed 
using SPH with sink particles replacing dense, collapsed cores. The different models are 
indicated in the figure, details can be found in table |I[ The figure shows how the efficiency 
of local collapse depends on the scale and strength of turbulent driving. 
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Fig. 5. — Mass fraction M,, in dense cores as function of time for the three 256^ ZEUS 
models driven with k = 1 — 2 (solid, crosses), A; = 3 — 4 (dashed, stars) and = 7 — 8 
(dotted, triangles). M* is the sum of all cores found by Clumpfind as discussed in the text. 
Note that the method identifies cores only after gravity is turned on, i.e. for t > 0.0. 
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Fig. 6. — Study of resolution and statistical variation of the core mass fraction M^, over time 
for SPH models with turbulent Jeans mass (Mj)turb = 3.2 and /c = 3 — 4. The low-resolution 
model B21 has 20 000 particles, for medium-resolution model B2 this number is 50 000, and 
for the high-resolution model B2h it is 200 000. Model B2 has been repeated four times with 
different realizations of the driving field. The alternative models B^"- to B^*^ are indicated 
by dotted fines. Note the large variance effect. 




Fig. 7. — Resolution study of core mass fraction as function of time for ZEUS models 
with turbulent Jeans mass (Mj)turb = 3.2 and driving wave number A; = 3 — 4. The models 
have resolutions of 64^ (dotted), 128"^ (dashed), and 256'^ cells (solid). is computed using 
Clumpfind as discussed in the text. 
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Fig. 8. — Energy as function of wave number k for models with different driving scale: (a) 
Al with k = 1 — 2, (b) A2 with k = 3 — 4 and (c) A3 with k = 7 — 8. The simulations 
are studied at t = 0.0, when the hydrodynamic turbulence is fully developed, immediately 
after gravity is included. The plots show potential energy Ep^t (thick solid lines), kinetic 
energy E'kin (thick long-dashed lines), its solenoidal component Esoi (short-dashed hnes) 
and its compressional component £'com (dotted lines). The thin dotted lines indicate the 
slope —5/3 expected from the Kolmogorov (1941) theory and the slope —2 expected for 
velocity discontinuities associated with shocks. For comparison, plot (d) shows a strongly 
self-gravitating model that completely lacks turbulent support and therefore contracts on all 
scales (data from Klessen et al. 1998). The energy spectra are computed on a 128^ grid onto 
which the SPH particle distribution has been assigned using the kernel smoothing procedure. 
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Fig. 9. — Fourier analysis of the high-resolution model B2h ((Mj)turb = 3.2 and A; = 3 — 4) 
at different stages of its dynamical evolution indicated on each plot. Notation and scaling 
are the same as in figure Again, the SPH model is sampled on a 128^ mesh. 
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Fig. 10. — Wave mode comparison between four models with identical physical parameters 
((Mj)turb = 3.2 and A; = 3 — 4) computed with different numerical methods and resolution: 
SPH models B2 and B2h with 50 000 and 200 000 particles, and ZEUS models V2 and V2h 
with 128^ and 256^ grid cells. To enable direct comparison, equivalent energy components 
of all four models are plotted in each panel. The upper half (a - d) of the figure shows the 
energy distribution of a state of fully developed hydrodynamic turbulence without gravity. 
The lower half (e - h) depicts the system after gravity is included, when M^. = 5% of the 
total mass is collapsed onto dense cores. Again, all spectra are computed on a grid with 
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Fig. 11. — Comparison of core locations between two globally stable models with different 
driving wave length (Bl with k = 1 — 2 and B3 with k = 7 — 8) projected into (a) the xy- 
plane and into (b) the a:2;-plane. Plots (c) and (d) show the core locations for model Bl now 
contrasted with a simulation of decaying turbulence from Klessen (1999). The snapshots 
arc selected such that the mass accumulated in dense cores is M* ^ 20%. Note the different 
times needed for the different models to reach this point. For model Bl data arc taken at 
t — 1.1, for S3 at i = 12.3. The simulation of freely decaying turbulence is shown at t — 1.1. 
All times are normahzed to the global free-fall time scale of the system. 
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Fig. 12. — Core positions for model Bl {k = 1 — 2) and the decay model when the core mass 



fraction is ^ 65%, projected into (a) the xy-plane and (b) the x2;-plane (cf. figure 11c 
& d). For Bl the time is t = 8.7 and for decay model t = 2.1. Whereas the cluster in 
Bl is completely dissolved and the stars are widely dispersed throughout the computational 
volume, the core cluster in the decay simulation remains bound. 
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Fig. 13. — Core masses as function of time in SPH models (a) Bl with k — 1 — 2 driving, 
(b) B2 with A; = 3 — 4 driving, and (c) BS with A; = 7 — 8 driving. The curves represent 
the formation and accretion histories of individual cores. For the sake of clarity, only every 
other core is shown in (a) and (b), whereas in (c) the evolution of every single core is plotted. 
Time is given in units of the global free-fall time rg. Note the different time scale in each 
plot. In the depicted time interval models Bl and B2 reach a core mass fraction M* = 70%, 
and both form roughly 50 cores. Model B3 reaches M^, = 35% and forms only 25 cores. 
Figure (d) compares the distributions of formation times. The age spread increases with 
decreasing driving scale showing that clustered core formation should lead to a coeval stellar 
population, whereas a distributed stellar population should exhibit considerable age spread. 



